import numpy as N
import pylab
plot = pylab.plot
import math

from scikits.learn.gmm import GMM, sample_gaussian

############################################################################## evaluate Gauss from Netlab

def gauss(x, sigma, mu):
	"""
	evaluate Gaussian density with mean MU and variance SIGMA
	not tested on 2 or higher dimension, only 1 dimension.
	
	based on Ian T Nabney (1996, 2001) 
	Netlab Toolkit for Matlab
	"""
	d = 1
	n = len(x)
	if sigma.ndim >=2:
		inverse_covariance = N.linalg.inv(sigma)
		determinant_covariance = N.linalg.det(sigma)
	else:
		inverse_covariance = [1.0/sigma]
		determinant_covariance = sigma
	x = x - N.ones(n)*mu
	fact = N.multiply((x*inverse_covariance),x)
	y = N.array([ math.exp(i) for i in fact[0]*-0.5])
	y = y/math.sqrt( (2*N.pi)**d * determinant_covariance  )
	return y	

############################################################################## Biased Gaussian Mixture Model
# TODO: test on this
class BiasedGMM(GMM):

	def __init__(self, *a, **k):
		print 'BiasedGMM __init__'
		GMM.__init__(self, *a, **k)

	def biasedSample(self, p, n):
		"""
		sample with probability p (0<p<1) from only n component
		"""

		#print 'BiasedSample > biasedSample'
		rand = N.random.rand()
		if rand < p:
			# sample biased
			if self._cvtype == 'tied':
				cv = self._covars
			else:
				cv = self._covars[n]
			return sample_gaussian(self._means[n], cv, self._cvtype)[0], 'biased' # [0] because sample_gaussian() return an array object
		else:
			# normal sample
			return self.rvs()[0][0], 'normal' # [0][0] because rsv() method return an array-shape object 

if __name__ == "__main__":

	# Assertation test.
	#print gauss(N.array([1,2,3,4,5]), N.array([[5]]), N.array([[4]]))
	# [ 0.07253707  0.11959342  0.16143423  0.17841241  0.16143423]

	from scikits.learn.gmm import GMM
	mu = N.array([[91.1855], [69.7003], [54.9596]])
	va = N.array([[262.1538], [50.3270], [9.5023]])
	prior = N.array([[0.3533], [0.3823], [0.2644]])
	gmm = GMM(n_dim=1, n_states=3,  means=mu, covars=va)
	
	plot(gmm.rvs(200))
	pylab.show()	
#	ymix = N.array([0 for i in range(150)])
#	x = N.array([j for j in range(150)])
#	for i in range(3):
#		y = gauss(x, va[i], mu[i])
#		ymix = ymix+y
#		plot(x, y, 'r')
#	plot(x, ymix, 'b')
#
#	pylab.show()
